library(knitr)
# purl("08_fec_multivariate.Rmd", output = "08_fec_multivariate.R", documentation = 2)

Setup

rerun_MCMC <- FALSE

set.seed(224819)

library(MCMCglmm)
## Loading required package: Matrix
## Loading required package: coda
## Loading required package: ape
library(tidyverse)
## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.0.0     ✔ purrr   0.2.5
## ✔ tibble  1.4.2     ✔ dplyr   0.7.6
## ✔ tidyr   0.8.1     ✔ stringr 1.3.1
## ✔ readr   1.1.1     ✔ forcats 0.3.0
## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(parallel)

source("color_map.R")
source("ggplot_theme.R")
h2fec <- read.table("../../Data/Processed/eggs_per_vial.txt",
                    sep = '\t',
                    dec = ".", header = TRUE,
                    stringsAsFactors = FALSE)

# Standardize egg_total
h2fec$egg_total <- as.numeric(scale(h2fec$egg_total))

h2fec$animal <- seq(1, nrow(h2fec))
h2fec$treat <- as.factor(h2fec$treat)

pedigree <- h2fec[, c("animal", "sireid", "damid")]
names(pedigree) <- c("animal", "sire", "dam")
pedigree$animal <- as.character(pedigree$animal)
pedigree$sire <- as.character(pedigree$sire)
pedigree$dam <- as.character(pedigree$dam)
sires <- data.frame(animal = unique(pedigree$sire),
                    sire = NA, dam = NA, stringsAsFactors = FALSE)
dams <- data.frame(animal = unique(pedigree$dam),
                   sire = NA, dam = NA, stringsAsFactors = FALSE)
pedigree <- bind_rows(sires, dams, pedigree) %>%
  as.data.frame()

genet_corr <- tibble(model = character(),
                     HS_STD = numeric(),
                     LY_STD = numeric(),
                     HS_LY = numeric(),
                     n_eff = numeric())

iter <- 6.5e6
burnin <- 5e4
thinning <- 500
chains <- 12

Multivariate analysis

if (rerun_MCMC) {
  HS <- h2fec %>% 
    filter(treat == "HS") %>% rename(Eggs_HS = egg_total) %>% 
    as.data.frame()
  DR <- h2fec %>% 
    filter(treat == "LY") %>% rename(Eggs_DR = egg_total) %>% 
    as.data.frame()
  STD <- h2fec %>%
    filter(treat == "STD") %>% rename(Eggs_STD = egg_total) %>% 
    as.data.frame()
  
  h2fec_mrg <- full_join(full_join(HS, DR), STD)
  
  prior1 <- list(R = list(V = diag(3) * 1.002, nu = 1.002),
                 G = list(G1 = list(V = diag(3) * 1.002, nu = 0.002)))
  
  priors <- list(prior1)
  prior_names <- c("Tri: V = diag(3) * 1.002, nu = 0.02")
  
  for (ii in 1:length(priors)) {
    prior <- priors[[ii]]
    fm <- mclapply(1:chains, function(i) {
      MCMCglmm(cbind(Eggs_STD, Eggs_HS, Eggs_DR) ~ trait - 1,
               random = ~ us(trait):animal,
               rcov = ~ idh(trait):units,
               data = h2fec_mrg,
               prior = prior,
               pedigree = pedigree,
               family = rep("gaussian", 3),
               nitt = iter,
               burnin = burnin,
               thin = thinning)
    }, mc.cores = chains)
    outfile <- paste0("../../Data/Processed/fec_total_multivariate_model_prior", ii, ".Rda")
    save(fm, file = outfile)
    
    re <- lapply(fm, function(m) m$VCV)
    re <- do.call(mcmc.list, re)
    re <- as.mcmc(as.matrix(re))
    
    n_eff <- effectiveSize(re)
    
    # STD vs. HS
    HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
      sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
             re[ , "traitEggs_HS:traitEggs_HS.animal"])
    
    # STD vs. DR
    DR_STD <- re[ , "traitEggs_DR:traitEggs_STD.animal"] /
      sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
             re[ , "traitEggs_DR:traitEggs_DR.animal"])
    
    # DR vs. HS
    HS_DR <- re[ , "traitEggs_HS:traitEggs_DR.animal"] /
      sqrt(re[ , "traitEggs_DR:traitEggs_DR.animal"] *
             re[ , "traitEggs_HS:traitEggs_HS.animal"])
    
    genet_corr <- bind_rows(genet_corr,
                            tibble(model = prior_names[[ii]],
                                   HS_STD = median(HS_STD),
                                   DR_STD = median(DR_STD),
                                   HS_DR = median(HS_DR),
                                   n_eff = mean(n_eff)))
  }
  
  write_csv(genet_corr, path = "../../Data/Processed/Genetic_Correlations_Fecundity.csv")
}

Analyze model

load("../../Data/Processed/fec_total_multivariate_model_prior1.Rda")

fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe[, 1, drop = FALSE], ask = FALSE)

plot(fe[, 2, drop = FALSE], ask = FALSE)

plot(fe[, 3, drop = FALSE], ask = FALSE)

# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)

plot(re[, 1, drop = FALSE], ask = FALSE)

plot(re[, 2, drop = FALSE], ask = FALSE)

plot(re[, 3, drop = FALSE], ask = FALSE)

autocorr.diag(re)
##           traitEggs_STD:traitEggs_STD.animal
## Lag 0                           1.0000000000
## Lag 500                         0.1874099388
## Lag 2500                        0.0368637247
## Lag 5000                        0.0068407312
## Lag 25000                      -0.0009086811
##           traitEggs_HS:traitEggs_STD.animal
## Lag 0                          1.0000000000
## Lag 500                        0.3135652226
## Lag 2500                       0.0488681915
## Lag 5000                       0.0004422247
## Lag 25000                      0.0066357983
##           traitEggs_DR:traitEggs_STD.animal
## Lag 0                           1.000000000
## Lag 500                         0.245970917
## Lag 2500                        0.044196514
## Lag 5000                        0.011865920
## Lag 25000                       0.002057849
##           traitEggs_STD:traitEggs_HS.animal
## Lag 0                          1.0000000000
## Lag 500                        0.3135652226
## Lag 2500                       0.0488681915
## Lag 5000                       0.0004422247
## Lag 25000                      0.0066357983
##           traitEggs_HS:traitEggs_HS.animal
## Lag 0                         1.0000000000
## Lag 500                       0.2983162660
## Lag 2500                      0.0387993500
## Lag 5000                      0.0007746715
## Lag 25000                    -0.0006467404
##           traitEggs_DR:traitEggs_HS.animal
## Lag 0                          1.000000000
## Lag 500                        0.329916759
## Lag 2500                       0.051896081
## Lag 5000                       0.002264740
## Lag 25000                      0.004759115
##           traitEggs_STD:traitEggs_DR.animal
## Lag 0                           1.000000000
## Lag 500                         0.245970917
## Lag 2500                        0.044196514
## Lag 5000                        0.011865920
## Lag 25000                       0.002057849
##           traitEggs_HS:traitEggs_DR.animal
## Lag 0                          1.000000000
## Lag 500                        0.329916759
## Lag 2500                       0.051896081
## Lag 5000                       0.002264740
## Lag 25000                      0.004759115
##           traitEggs_DR:traitEggs_DR.animal traitEggs_STD.units
## Lag 0                          1.000000000        1.0000000000
## Lag 500                        0.343334829        0.1694615436
## Lag 2500                       0.071567781        0.0343858227
## Lag 5000                       0.012262294        0.0066624941
## Lag 25000                     -0.003595793       -0.0006694499
##           traitEggs_HS.units traitEggs_DR.units
## Lag 0           1.0000000000       1.0000000000
## Lag 500         0.1354468964       0.1653264427
## Lag 2500        0.0181847424       0.0320639066
## Lag 5000       -0.0004421368       0.0057609407
## Lag 25000      -0.0045816742       0.0007743094
effectiveSize(re)
## traitEggs_STD:traitEggs_STD.animal  traitEggs_HS:traitEggs_STD.animal 
##                           74878.16                           59041.49 
##  traitEggs_DR:traitEggs_STD.animal  traitEggs_STD:traitEggs_HS.animal 
##                           64199.26                           59041.49 
##   traitEggs_HS:traitEggs_HS.animal   traitEggs_DR:traitEggs_HS.animal 
##                           64269.67                           57406.62 
##  traitEggs_STD:traitEggs_DR.animal   traitEggs_HS:traitEggs_DR.animal 
##                           64199.26                           57406.62 
##   traitEggs_DR:traitEggs_DR.animal                traitEggs_STD.units 
##                           50166.62                           78571.11 
##                 traitEggs_HS.units                 traitEggs_DR.units 
##                           96791.83                           78128.78
# Concatenate samples
re <- as.mcmc(as.matrix(re))

head(re)
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1 
## End = 7 
## Thinning interval = 1 
##      traitEggs_STD:traitEggs_STD.animal traitEggs_HS:traitEggs_STD.animal
## [1,]                          0.2730985                        0.12469556
## [2,]                          0.2995462                        0.07412571
## [3,]                          0.4452922                        0.18096616
## [4,]                          0.1252125                        0.07514900
## [5,]                          0.1048659                        0.04575423
## [6,]                          0.1438339                        0.05751000
## [7,]                          0.2007193                        0.03386499
##      traitEggs_DR:traitEggs_STD.animal traitEggs_STD:traitEggs_HS.animal
## [1,]                        0.24849777                        0.12469556
## [2,]                        0.20893779                        0.07412571
## [3,]                        0.20795884                        0.18096616
## [4,]                        0.11505847                        0.07514900
## [5,]                        0.08302472                        0.04575423
## [6,]                        0.08359697                        0.05751000
## [7,]                        0.13166049                        0.03386499
##      traitEggs_HS:traitEggs_HS.animal traitEggs_DR:traitEggs_HS.animal
## [1,]                      0.110130511                       0.11774301
## [2,]                      0.034842647                       0.04079747
## [3,]                      0.118196674                       0.07136862
## [4,]                      0.070105046                       0.08114866
## [5,]                      0.022681420                       0.03957909
## [6,]                      0.035602734                       0.03912530
## [7,]                      0.008878417                       0.02152333
##      traitEggs_STD:traitEggs_DR.animal traitEggs_HS:traitEggs_DR.animal
## [1,]                        0.24849777                       0.11774301
## [2,]                        0.20893779                       0.04079747
## [3,]                        0.20795884                       0.07136862
## [4,]                        0.11505847                       0.08114866
## [5,]                        0.08302472                       0.03957909
## [6,]                        0.08359697                       0.03912530
## [7,]                        0.13166049                       0.02152333
##      traitEggs_DR:traitEggs_DR.animal traitEggs_STD.units
## [1,]                       0.23588453           0.4891190
## [2,]                       0.16338372           0.4610870
## [3,]                       0.10745186           0.4486223
## [4,]                       0.11301484           0.6035127
## [5,]                       0.07074388           0.6815414
## [6,]                       0.05895181           0.4694439
## [7,]                       0.09380421           0.5364391
##      traitEggs_HS.units traitEggs_DR.units
## [1,]          0.1884776          0.2815059
## [2,]          0.1863606          0.2473274
## [3,]          0.1944373          0.3105254
## [4,]          0.1768785          0.3693298
## [5,]          0.2189036          0.3670518
## [6,]          0.2347433          0.4200357
## [7,]          0.2539427          0.3827891
# STD vs. HS
plot(re[ , "traitEggs_HS:traitEggs_STD.animal"])

plot(re[ , "traitEggs_STD:traitEggs_STD.animal"])

HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
  sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
         re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_STD)

median(HS_STD)
## [1] 0.7124564
HPDinterval(HS_STD)
##           lower     upper
## var1 -0.5112175 0.9987634
## attr(,"Probability")
## [1] 0.95
# STD vs. DR
DR_STD <- re[ , "traitEggs_DR:traitEggs_STD.animal"] /
  sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
         re[ , "traitEggs_DR:traitEggs_DR.animal"])
plot(DR_STD)

median(DR_STD)
## [1] 0.9312949
HPDinterval(DR_STD)
##          lower   upper
## var1 0.3639826 0.99934
## attr(,"Probability")
## [1] 0.95
# DR vs. HS
HS_DR <- re[ , "traitEggs_HS:traitEggs_DR.animal"] /
  sqrt(re[ , "traitEggs_DR:traitEggs_DR.animal"] *
         re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_DR)

median(HS_DR)
## [1] 0.7238309
HPDinterval(HS_DR)
##           lower     upper
## var1 -0.5517535 0.9979309
## attr(,"Probability")
## [1] 0.95
save(re, file = "../../Data/Processed/fec_total_multivariate_model_output.Rda")

Plot for paper

library(tidyverse)
library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
B <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
                `DR vs. STD` = as.numeric(DR_STD),
                `HS vs. DR` = as.numeric(HS_DR))

B <- as_tibble( B ) %>%  
  select(`HS vs. STD`, `DR vs. STD`, `HS vs. DR`) %>% 
  rename(`HS vs. C` = `HS vs. STD`) %>% 
  rename(`DR vs. C` = `DR vs. STD`) %>% 
  rename(`HS vs. DR` = `HS vs. DR`)

colMeans(B)
##  HS vs. C  DR vs. C HS vs. DR 
## 0.5413373 0.8384703 0.5400363
B %>% gather(Comparison, value) %>%
  ggplot(aes(value, color = Comparison)) +
  geom_line(stat = "density", size = 1) +
  labs(x = "Genetic Correlation", y = "Density") +
  theme(legend.position = c(0.10, 0.85),
        text = element_text(size = 10),
        legend.text = element_text(size = 10))+
  scale_x_continuous(expand = c(0, 0), limits = c(-1.2, 1.2)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 12)) +
  my_theme

ggsave(last_plot(), file = "../../Figures/Fecundity_Genetic_Correlation_Plot.pdf",
       width = 4, height = 4)
Fecundity_correlation <- last_plot()
save(Fecundity_correlation, file = "../../Figures/Fecundity_correlation.Rda")

Parameter expanded multivariate analysis

if (rerun_MCMC) {
  
  iter <- 1e5
  burnin <- 2e3
  thinning <- 500
  chains <- 12

  HS <- h2fec %>% 
    filter(treat == "HS") %>% rename(Eggs_HS = egg_total) %>% 
    as.data.frame()
  DR <- h2fec %>% 
    filter(treat == "LY") %>% rename(Eggs_DR = egg_total) %>% 
    as.data.frame()
  STD <- h2fec %>%
    filter(treat == "STD") %>% rename(Eggs_STD = egg_total) %>% 
    as.data.frame()
  
  h2fec_mrg <- full_join(full_join(HS, DR), STD)
  
  prior1 <- list(R = list(V = diag(3), nu = 1.002, fix = 3),
                 G = list(G1 = list(V = diag(3),
                                    nu = 1,
                                    alpha.mu = rep(0, 3),
                                    alpha.V = diag(3) * 25)))
  
  priors <- list(prior1)
  prior_names <- c("Tri: Parameter expanded")
  
  for (ii in 1:length(priors)) {
    prior <- priors[[ii]]
    fm <- mclapply(1:chains, function(i) {
      MCMCglmm(cbind(Eggs_STD, Eggs_HS, Eggs_DR) ~ trait - 1,
               random = ~ us(trait):animal,
               rcov = ~ idh(trait):units,
               data = h2fec_mrg,
               prior = prior,
               pedigree = pedigree,
               family = rep("gaussian", 3),
               nitt = iter,
               burnin = burnin,
               thin = thinning)
    }, mc.cores = chains)
    outfile <- paste0("../../Data/Processed/fec_total_multivariate_model_prior_expanded", ii, ".Rda")
    save(fm, file = outfile)
    
    re <- lapply(fm, function(m) m$VCV)
    re <- do.call(mcmc.list, re)
    re <- as.mcmc(as.matrix(re))
    
    n_eff <- effectiveSize(re)
    
    # STD vs. HS
    HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
      sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
             re[ , "traitEggs_HS:traitEggs_HS.animal"])
    
    # STD vs. DR
    DR_STD <- re[ , "traitEggs_DR:traitEggs_STD.animal"] /
      sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
             re[ , "traitEggs_DR:traitEggs_DR.animal"])
    
    # DR vs. HS
    HS_DR <- re[ , "traitEggs_HS:traitEggs_DR.animal"] /
      sqrt(re[ , "traitEggs_DR:traitEggs_DR.animal"] *
             re[ , "traitEggs_HS:traitEggs_HS.animal"])
    
    genet_corr <- bind_rows(genet_corr,
                            tibble(model = prior_names[[ii]],
                                   HS_STD = median(HS_STD),
                                   DR_STD = median(DR_STD),
                                   HS_DR = median(HS_DR),
                                   n_eff = mean(n_eff)))
  }
  
  write_csv(
    genet_corr,
    path = "../../Data/Processed/Genetic_Correlations_Fecundity_Parameter_Expanded.csv")
}

Analyze parameter expanded model

load("../../Data/Processed/fec_total_multivariate_model_prior_expanded1.Rda")

fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe[, 1, drop = FALSE], ask = FALSE)

plot(fe[, 2, drop = FALSE], ask = FALSE)

plot(fe[, 3, drop = FALSE], ask = FALSE)

# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)

plot(re[, 1, drop = FALSE], ask = FALSE)

plot(re[, 2, drop = FALSE], ask = FALSE)

plot(re[, 3, drop = FALSE], ask = FALSE)

autocorr.diag(re)
##           traitEggs_STD:traitEggs_STD.animal
## Lag 0                            1.000000000
## Lag 500                         -0.003725192
## Lag 2500                         0.019490364
## Lag 5000                         0.030377401
## Lag 25000                       -0.001895760
##           traitEggs_HS:traitEggs_STD.animal
## Lag 0                           1.000000000
## Lag 500                        -0.016525153
## Lag 2500                       -0.013641051
## Lag 5000                        0.002471646
## Lag 25000                       0.010064237
##           traitEggs_DR:traitEggs_STD.animal
## Lag 0                           1.000000000
## Lag 500                         0.028389065
## Lag 2500                        0.006133787
## Lag 5000                       -0.040036034
## Lag 25000                      -0.001075109
##           traitEggs_STD:traitEggs_HS.animal
## Lag 0                           1.000000000
## Lag 500                        -0.016525153
## Lag 2500                       -0.013641051
## Lag 5000                        0.002471646
## Lag 25000                       0.010064237
##           traitEggs_HS:traitEggs_HS.animal
## Lag 0                         1.0000000000
## Lag 500                       0.0403265957
## Lag 2500                      0.0024874149
## Lag 5000                     -0.0045582535
## Lag 25000                     0.0003642594
##           traitEggs_DR:traitEggs_HS.animal
## Lag 0                         1.0000000000
## Lag 500                      -0.0108965197
## Lag 2500                      0.0198383220
## Lag 5000                     -0.0001302109
## Lag 25000                     0.0080958169
##           traitEggs_STD:traitEggs_DR.animal
## Lag 0                           1.000000000
## Lag 500                         0.028389065
## Lag 2500                        0.006133787
## Lag 5000                       -0.040036034
## Lag 25000                      -0.001075109
##           traitEggs_HS:traitEggs_DR.animal
## Lag 0                         1.0000000000
## Lag 500                      -0.0108965197
## Lag 2500                      0.0198383220
## Lag 5000                     -0.0001302109
## Lag 25000                     0.0080958169
##           traitEggs_DR:traitEggs_DR.animal traitEggs_STD.units
## Lag 0                          1.000000000         1.000000000
## Lag 500                        0.023406187        -0.022469930
## Lag 2500                       0.010900801         0.013481221
## Lag 5000                      -0.004557295         0.020502363
## Lag 25000                      0.011769764        -0.006948712
##           traitEggs_HS.units traitEggs_DR.units
## Lag 0            1.000000000                NaN
## Lag 500         -0.005884207                NaN
## Lag 2500         0.008485809                NaN
## Lag 5000         0.014211774                NaN
## Lag 25000       -0.010895854                NaN
effectiveSize(re)
## traitEggs_STD:traitEggs_STD.animal  traitEggs_HS:traitEggs_STD.animal 
##                           2437.161                           2464.272 
##  traitEggs_DR:traitEggs_STD.animal  traitEggs_STD:traitEggs_HS.animal 
##                           2381.271                           2464.272 
##   traitEggs_HS:traitEggs_HS.animal   traitEggs_DR:traitEggs_HS.animal 
##                           2444.515                           2431.644 
##  traitEggs_STD:traitEggs_DR.animal   traitEggs_HS:traitEggs_DR.animal 
##                           2381.271                           2431.644 
##   traitEggs_DR:traitEggs_DR.animal                traitEggs_STD.units 
##                           2303.084                           2460.700 
##                 traitEggs_HS.units                 traitEggs_DR.units 
##                           2372.235                              0.000
# Concatenate samples
re <- as.mcmc(as.matrix(re))

head(re)
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1 
## End = 7 
## Thinning interval = 1 
##      traitEggs_STD:traitEggs_STD.animal traitEggs_HS:traitEggs_STD.animal
## [1,]                          0.3590307                      0.1282773450
## [2,]                          0.2442366                      0.0006156665
## [3,]                          0.1861954                     -0.0100121837
## [4,]                          0.4102305                      0.1024781219
## [5,]                          0.1929804                      0.0546066988
## [6,]                          0.2738059                      0.0193737220
## [7,]                          0.2613179                      0.0017872135
##      traitEggs_DR:traitEggs_STD.animal traitEggs_STD:traitEggs_HS.animal
## [1,]                     -0.0083137086                      0.1282773450
## [2,]                      0.1011918879                      0.0006156665
## [3,]                      0.0018253463                     -0.0100121837
## [4,]                      0.0001729217                      0.1024781219
## [5,]                      0.0564094448                      0.0546066988
## [6,]                     -0.0117165610                      0.0193737220
## [7,]                      0.0178162082                      0.0017872135
##      traitEggs_HS:traitEggs_HS.animal traitEggs_DR:traitEggs_HS.animal
## [1,]                     7.439954e-02                    -0.0051252069
## [2,]                     1.955109e-05                     0.0005158809
## [3,]                     6.917970e-02                    -0.0003176388
## [4,]                     5.379830e-02                     0.0036724221
## [5,]                     8.522121e-02                     0.0503100316
## [6,]                     4.261395e-02                     0.0069020223
## [7,]                     6.829332e-02                     0.0013074032
##      traitEggs_STD:traitEggs_DR.animal traitEggs_HS:traitEggs_DR.animal
## [1,]                     -0.0083137086                    -0.0051252069
## [2,]                      0.1011918879                     0.0005158809
## [3,]                      0.0018253463                    -0.0003176388
## [4,]                      0.0001729217                     0.0036724221
## [5,]                      0.0564094448                     0.0503100316
## [6,]                     -0.0117165610                     0.0069020223
## [7,]                      0.0178162082                     0.0013074032
##      traitEggs_DR:traitEggs_DR.animal traitEggs_STD.units
## [1,]                      0.004483110           0.4588685
## [2,]                      0.058540298           0.3533715
## [3,]                      0.000445515           0.5218645
## [4,]                      0.001503260           0.1974659
## [5,]                      0.046957532           0.6654494
## [6,]                      0.003940158           0.4780136
## [7,]                      0.001747412           0.4619696
##      traitEggs_HS.units traitEggs_DR.units
## [1,]          0.2339418                  1
## [2,]          0.3095291                  1
## [3,]          0.1936399                  1
## [4,]          0.2018433                  1
## [5,]          0.1806931                  1
## [6,]          0.1818133                  1
## [7,]          0.1943753                  1
# STD vs. HS
plot(re[ , "traitEggs_HS:traitEggs_STD.animal"])

plot(re[ , "traitEggs_STD:traitEggs_STD.animal"])

HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
  sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
         re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_STD)

median(HS_STD)
## [1] 0.3946399
HPDinterval(HS_STD)
##           lower     upper
## var1 -0.7472588 0.9999947
## attr(,"Probability")
## [1] 0.9498299
# STD vs. DR
DR_STD <- re[ , "traitEggs_DR:traitEggs_STD.animal"] /
  sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
         re[ , "traitEggs_DR:traitEggs_DR.animal"])
plot(DR_STD)

median(DR_STD)
## [1] 0.405907
HPDinterval(DR_STD)
##           lower     upper
## var1 -0.8619617 0.9999912
## attr(,"Probability")
## [1] 0.9498299
# DR vs. HS
HS_DR <- re[ , "traitEggs_HS:traitEggs_DR.animal"] /
  sqrt(re[ , "traitEggs_DR:traitEggs_DR.animal"] *
         re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_DR)

median(HS_DR)
## [1] 0.2382352
HPDinterval(HS_DR)
##           lower     upper
## var1 -0.9068489 0.9999722
## attr(,"Probability")
## [1] 0.9498299

Parameter expanded plot

library(tidyverse)
library(cowplot)

B <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
                `DR vs. STD` = as.numeric(DR_STD),
                `HS vs. DR` = as.numeric(HS_DR))

B <- as_tibble( B ) %>%  
  select(`HS vs. STD`, `DR vs. STD`, `HS vs. DR`) %>% 
  rename(`HS vs. C` = `HS vs. STD`) %>% 
  rename(`DR vs. C` = `DR vs. STD`) %>% 
  rename(`HS vs. DR` = `HS vs. DR`)

colMeans(B)
##  HS vs. C  DR vs. C HS vs. DR 
## 0.2945234 0.2643091 0.1457471
B %>% gather(Comparison, value) %>%
  ggplot(aes(value, color = Comparison)) +
  geom_line(stat = "density", size = 1) +
  labs(x = "Genetic Correlation", y = "Density") +
  theme(legend.position = c(0.10, 0.85),
        text = element_text(size = 10),
        legend.text = element_text(size = 10))+
  scale_x_continuous(expand = c(0, 0), limits = c(-1, 1)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 12)) +
  my_theme

ggsave(last_plot(), file = "../../Figures/Fecundity_Genetic_Correlation_Plot.pdf",
       width = 4, height = 4)
Fecundity_correlation <- last_plot()
save(Fecundity_correlation, file = "../../Figures/Fecundity_correlation.Rda")